﻿using System;
using System.Collections.Generic;
using GeoAPI.CoordinateSystems;
using GeoAPI.CoordinateSystems.Transformations;

namespace ProjNet.CoordinateSystems.Projections
{
    public class PolarStereoNorth : MapProjection
    {
        private const double DEG2RAD = 0.017453292519943296;
        private const double RAD2DEG = 57.295779513082322;
        private const double PI = 3.14159265358979323846;
        private const double A = 6378137; //WGS84 // +a=6378137 +b=6356752.3142
        private const double E = 0.0818191913108718138;
        private const double K = 1;
        private double _false_easting = 0;
        private double _false_northing = 0;

        public PolarStereoNorth(List<ProjectionParameter> parameters)
            : this(parameters, false)
        {
        }

        public PolarStereoNorth(List<ProjectionParameter> parameters, bool isInverse)
            : base(parameters, isInverse)
        {
            Name = "Polar Stereographic North";
            Authority = "AzukiMap";
            AuthorityCode = 9810;
            ProjectionParameter falseEasting = GetParameter("false_easting");
            ProjectionParameter falseNorthing = GetParameter("false_northing");

            if (falseEasting == null)
                throw new ArgumentException("Missing projection parameter 'false_easting'");
            if (falseNorthing == null)
                throw new ArgumentException("Missing projection parameter 'false_northing'");

            _false_easting = falseEasting.Value;
            _false_northing = falseNorthing.Value;
        }

        public override double[] MetersToDegrees(double[] p)
        {
            double e = p[0];
            double n = p[1];
            double lo;
            if (e > 0 && n > 0)
            {
                lo = 180;
            }
            else if (e < 0 && n > 0)
            {
                lo = -180;
            }
            else
            {
                lo = 0;
            }

            double lon = 0 * DEG2RAD + Math.Atan(((e - _false_easting) / (_false_northing - n)));

            double x = Math.PI / 2 - 2 * Math.Atan(Tt(e, n));
            double lat = x +
                         ((Math.Pow(E, 2) / 2) + (5 * Math.Pow(E, 4) / 24) + (Math.Pow(E, 6) / 12) + (13 * Math.Pow(E, 8) / 360)) *
                         Math.Sin(2 * x) +
                         (7 * Math.Pow(E, 4) / 48 + 29 * Math.Pow(E, 6) / 240 + 811 * Math.Pow(E, 8) / 11520) * Math.Sin(4 * x) +
                         (7 * Math.Pow(E, 6) / 120 + 81 * Math.Pow(E, 8) / 1120) * Math.Sin(6 * x) +
                         (4279 * Math.Pow(E, 8) / 161280) * Math.Sin(8 * x);
            lon = lo + (RAD2DEG * lon);
            lat = RAD2DEG * lat;
            lat = (lat == 90) ? 89.999999 : lat;
            return new[] { lon, lat };
        }

        public override double[] DegreesToMeters(double[] lonlat)
        {
            if (double.IsNaN(lonlat[0]) || double.IsNaN(lonlat[1]))
                return new[] { Double.NaN, Double.NaN, };

            double dLon = Degrees2Radians(AdjustLongitude(lonlat[0]));
            double dLat = Degrees2Radians(lonlat[1]);

            const double lo = 0;
            double e = _false_easting + Rho(dLat) * Math.Sin(dLon - lo);
            double n = _false_northing - Rho(dLat) * Math.Cos(dLon - lo);

            return new[] { e, n };
        }

        public override IMathTransform Inverse()
        {
			if (_inverse == null)
				_inverse = new PolarStereoNorth(this._Parameters, ! _isInverse);			
			return _inverse;
        }

        private static double Rho(double lat)
        {
            return 2 * A * K * T(lat) / Math.Pow(Math.Pow((1 + E), (1 + E)) * Math.Pow((1 - E), (1 - E)), 0.5);
        }

        private static double T(double lat)
        {
            return Math.Tan((PI / 4) - (lat / 2)) / Math.Pow(((1 - (E * Math.Sin(lat))) / (1 + (E * Math.Sin(lat)))), E / 2);
        }

        private double Tt(double e, double n)
        {
            double p = Pp(e, n);
            return p * Math.Pow((Math.Pow(1 + E, 1 + E) * Math.Pow(1 - E, 1 - E)), 0.5) / (2 * A * K);
        }

        private double Pp(double e, double n)
        {
            return Math.Sqrt(Math.Pow((e - _false_easting), 2) + Math.Pow((n - _false_northing), 2));
        }

        private double AdjustLongitude(double lon )
        {
            if (-180 <= lon && lon <= 180)
                return lon;
            if ( lon < -180 )
            {
                lon += 180;
                return AdjustLongitude(lon);
            }
            if ( 180 < lon)
            {
                lon -= 180;
                return AdjustLongitude(lon);
            }
            // never come here.
            return 0;
        }
    }
}